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Abstract 



/^ ', A typical procedure to integrate numerically the time dependent Schro- 

^H ■ dinger equation involves two stages. In the first one carries out a space 

(-H , discretization of the continuous problem. This results in the linear system 

of differential equations idu/dt = Hu, where H is a real symmetric matrix, 
whose solution with initial value u(0) = mq G C^ is given by u{t) = 
e~^*^uo. Usually, this exponential matrix is expensive to evaluate, so that 
time stepping methods to construct approximations to u from time t„ to 
O^l ' tn+i arc considered in the second phase of the procedure. Among them, 

^ I schemes involving multiplications of the matrix H with vectors, such as 

s!j ' Lanczos and Chebyshcv methods, are particularly efficient. 

r — , In this work we consider a particular class of splitting methods which 

■^^ ' also involves only products Hu. We carry out an error analysis of these 

integrators and propose a strategy which allows us to construct difFer- 
f~^ ' ent splitting symplectic methods of different order (even of order zero) 

f^ I possessing a large stability interval that can be adapted to different space 

regularity conditions and different accuracy ranges of the spatial discretiza- 
tion. The validity of the procedure and the performance of the resulting 
schemes are illustrated on several numerical examples. 
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1 Introduction 

To describe and understand the dynamics and evolution of many basic atomic 
and molecular phenomena, their time dependent quantum mechanical treat- 
ment is essential. Thus, for instance, in molecular dynamics, the construction 
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of models and simulations of molecular encounters can benefit a good deal from 
time dependent computations. The same applies to scattering processes such 
as atom-diatom collisions and triatomic photo-dissociation and, in general, to 
quantum mechanical phenomena where there is an initial state that under the 
influence of a given potential evolves through time to achieve a final asymptotic 
state (e.g., chemical reactions, unimolecular breakdown, desorption, etc.). This 
requires, of course, to solve the time dependent Schrodinger equation (h = 1) 

i-tP{x,t) = Hijj{x,t), (1) 

where H is the Hamiltonian operator, ip : R x M — > C is the wave function 
representing the state of the system and the initial state is ip{x,0) = 'ipo{x). 
Usually 

H = f{P) + ViX) = ^P2 ^ y^x-^ (2) 

and the operators X, P are defined by their actions on ip{x,t) as 

Xtp{x,t) = xip{x,t), P ip{x,t) = -iV^{x,t). 

The solution of ([T]) provides all dynamical information on the physical system 
at any time. It can be expressed as 

V'(x,t) = [/(t)Vo(x), (3) 

where U represents the evolution operator, which is linear and satisfies the 
equation idU{t)/dt = HU{t) with C/(0) = /. Since the Hamiltonian is explicitly 
time independent, the evolution operator is given formally by 

lJ{t) = e"^*^. (4) 

In practice, however, the Schrodinger equation has to be solved numerically, 
and the procedure involves basically two steps. The first one consists in consid- 
ering a faithful discrete spatial representation of the initial wave function '4^o{x) 
and the operator H on an appropriately constructed grid. Once this spatial 
discretization is built, the initial wave function is propagated in time until the 
end of the dynamical event. It is on the second stage of this process where we 
will concentrate our analysis. 

As a result of the discretization of eq. ([TJ in space, one is left with a lin- 
ear equation idu/dt = Hu, with a Hermitian matrix H of large dimension 
and large norm. Since evaluating exactly the exponential exp{—itH) is com- 
putationally expensive, approximation methods requiring only matrix-vector 
products with H are particularly appropriate [22]. Among them, the class of 
splitting symplectic methods has received considerable attention in the litera- 
ture [131 [231 [271 [21 [3] . Ill this case exp{—itH) is approximated by a composition 
of symplectic matrices. While it has been shown that stable high order methods 
belonging to this family do exist, such high degree of accuracy may be dispro- 
portionate in comparison with the error involved in the spatial discretization. 



and also inappropriate particularly when the problem at hand involves non- 
smooth solutions, as high order methods make small phase errors in the low 
frequencies but much larger errors in the high frequencies. An error analysis of 
this family of integrators, in particular, could help one to design different effi- 
cient time integrators adapted to different accuracy requirements and spacial 
regularity situations. 

The analysis carried out in the present paper could be considered a step 
forward in this direction. We present a strategy which allows us to construct 
different splitting symplectic methods of different order and large stability in- 
terval (with a large number of stages) that can be adapted to different space 
regularity conditions and different accuracy ranges of the spatial discretization. 
When this regularity degree is low, sometimes the best option is provided by 
methods of order zero. 

Since the splitting methods we analyze here only involve products of the 
matrix H with vectors, they belong to the same class of integrators as the 
Chebyshev and Lanczos methods, in the sense that all of them approximate 
exp{—itH)uQ by linear combinations of terms of the form H^uq (j > 1). 

The plan of the paper is as follows. In section [2] we review first the Fourier 
collocation approach carrying out the spatial discretization of the Schrodinger 
equation, and then we turn our attention to the time discretization errors of 
symplectic splitting methods. The bulk of the paper is contained in section 
[31 There we carry out a theoretical analysis of symplectic splitting methods 
and obtain some estimates on the time discretization error. These estimates 
in turn allows us to build different classes of splitting schemes in section U 
which are then illustrated in section [5] on several numerical examples exhibiting 
different degrees of regularity. Here we also include, for comparison, results 
achieved by the Lanczos and Chebyshev methods. Finally, section [6] contains 
some conclusions. 

2 Space and time discretization 

Among many possible ways to discretize the Schrodinger equation in space, 
collocation spectral methods possess several attractive features: they allow a 
relatively small grid size for representing the wave function, are simple to im- 
plement and provide an extremely high order of accuracy if the solution of the 
problem is sufficiently smooth [ll^ I12j. In fact, spectral methods are superior 
to local methods (such as finite difference schemes) not only when very high 
spatial resolution is required, but also when long time integration is carried out, 
since the resulting spatial discretization does not cause a deterioration of the 
phase error as the integration in time goes on |16j . 

To simplify the treatment, we will limit ourselves to the one-dimensional 
case and assume that the wave function is negligible outside an interval [a,/?]. 
In such a situation one may reformulate the problem on the finite interval with 
periodic boundary conditions. After rescaling, one may assume without loss of 



i— ^(x,t) = - — Tr^(x,t) + T/(x)V(x,t), 0<x<27r (5) 



generality that the space interval is [0,27r], and therefore 

with ■ilj{0,t) = i/;(27r,t) for all t. In the Fourier-collocation (or pseudospectral) 

approach, one intends to construct approximations based on the equidistant 

interpolation grid 

27r 
Xj = j^j, j = 0,...,N -1 

where A'^ is even (although the formalism can also be adapted to an odd number 
of points). Then one seeks a solution of the form j22] 



V'^(x,t)= Yl Cn(t)e*"", xG[0,27r) (6) 

|n|<Ar/2 

where the coefficients Cn{t) are related to the grid values ijJN{xj,t) through 
a discrete Fourier transform of length A'^, J-'n [26]. Its computation can be 
accomplished by the Fast Fourier Transform (FFT) algorithm with 0{N log N) 
floating point operations. 

In the collocation approach, the grid values i{jN{xj,t) are determined by 
requiring that the approximation ([6]) satisfies the Schrodinger equation precisely 
at the grid points Xj |22j . This yields a system of A^ ordinary differential 
equations to determine the N point values 'ipi\f{xj,t): 



dt 
where 



i^- = Tj^^DnTnu + Vnu = Hu, u = {uo,ui,...,UN~i), (7) 



1 - -.2, 



Dn = — diag(n^), Vn = diag(y(xj)) (8) 

for n = —N/2, . . . , N/2 — 1 and j = 0, . . . , iV — 1. Observe that the matrices on 
the right-hand side of ([7|) are Hermitian. 

An important qualitative feature of this space discretization procedure is 
that it replaces the original Hilbert space /I^(0,27r) defined by the quantum 
mechanical problem by a discrete one in which the action of operators are 
approximated by A^ x A^ (Hermitian) matrices obeying the same quantum me- 
chanical commutation relations [18]. From a quantitative point of view, if the 
function tp is sufficiently smooth and periodic, then the coefficients Cn exhibit 
a rapid decay (in some cases, faster than algebraically in n~^, uniformly in A^), 
so that typically the value of A^ in the expansion ([6]) needs not to be very large 
to represent accurately the solution. Specifically, in [22] the following result is 
proved. 

Theorem 1 Suppose that the exact solution ■ip{x,t) of ^ is such that, for 
some s > 1, d^'^'^'ip{-,t) G 11^(0, 2tt) for every t > 0. Then the error due to the 
approximation 'il)i\f{x,t) defined by 1^ in the collocation approach is hounded by 

||V';v(-,t)-^(-,t)ll<CA^-^(l + t) maxJ|9^+V(-,t')||, 

where C depends only on s. 



When the problem is not periodic, the use of a truncated Fourier series 
introduces errors in the computation. In that case several techniques have been 
proposed to minimize its effects (see [TJ [S] and references therein) . 

The previous treatment can be generalized to several spatial dimensions, 
still exploiting all the one-dimensional features, by taking tensor products of 
one-dimensional expansions. The resulting functions are then defined on the 
Cartesian product of intervals [H [22] . 

We can then conclude that after the previous space discretization has been 
applied to eq. ([5]), one ends up with a linear system of ODEs of the form 

i—u{t) = Hu{t), n(0) = no G C^, (9) 

where H is a real symmetric matrix. This is the starting point for carrying out 
an integration in time. Although a collocation approach has been applied here, 
in fact any space discretization scheme leading to an equation of the form ([9]) 
fits in our subsequent analysis. 

The spatial discretization chosen has of course a direct consequence on the 
time propagation of the (discrete) wave function u{t), since the matrix H repre- 
senting the Hamiltonian has a discrete spectrum which depends on the scheme. 
This discrete representation, in addition, restricts the energy range of the prob- 
lem and therefore imposes an upper bound to the high frequency components 
represented in the propagation |19] . 

The exact solution of eq. ^ is given by 

u{t) = e-''^" uq, (10) 

but to compute the exponential of the N x N complex and full matrix —itH 
(typically also of large norm) by diagonalizing the matrix H can be prohibitively 
expensive for large values of A^. In practice, thus, one turns to time stepping 
methods advancing the approximate solution from time t„ to tn+i = t„ , + Ai, so 
that the aim is to construct an approximation Un+i ~ u{tn+i) = e~^ u{tn) 
as a map Un+i = 4>AtUn- 

Among them, exponential splitting schemes have been widely used when 
the Hamiltonian has the form given by ([2]) [9l [T9| I22j . In that case, equation 
dl]) reads 

iu={T + V)u, u{Q) = uo, (11) 

where ^ is a diagonal matrix associated with V and T is related to the kinetic 
energy T. It turns out that the solutions e~**^tio and e~**^no of equations 
iii = Tu and iu = Vu, respectively, can be easily determined [22], so that one 
may consider compositions of the form 

-ibmrV -iamTT _ _ _ -ibirV -iairT /-^2) 

where r = At. In (|12p the number of exponentials m (and therefore the number 
of coefficients {ai,bi}^i) has to be sufficiently large to solve all the equations 
required to achieve order r (the so called order conditions). 

Splitting methods of this class have several structure-preserving properties. 
They are unitary, so that the norm of u is preserved along the integration, and 



time-reversible when the composition (J12p is symmetric. Error estimates of 
such methods apphed to the Schrodinger equation [T71 [2ll [25] seem to suggest 
that, while they are indeed very efficient for high spatial regularity, they may 
not be very appropriate under conditions of limited regularity. 

Here we will concentrate on another class of splitting methods that have 
been considered in the literature [131 ESI HZl [21 [3] . Notice that the corresponding 
H in eq. ([7]) is a real symmetric matrix, and thus e~ is not only unitary, but 
also symplectic with canonical coordinates q = Re(u) and momenta p = Im(u). 
In consequence, equation ([9]) is equivalent to [HI [H] 

q = Hp, p = —Hq. (13) 

Alternatively, one may write 

i(:)K-". o)(:)^(--)(:)> <") 

with the 2N x 2N matrices A and B given by 

^='o ' ^= -F 



The solution operator corresponding to (J14p can be written in terms of the 
rotation matrix 

^^y)-[-sm{y) cos{y) ) ^^^> 

as 0{tH), which is an orthogonal and symplectic 2N x 2A^ matrix. Computing 
0{tH) exactly (by diagonalizing the matrix H) is, as mentioned before for 
its complex representation e~**^, computationally very expensive, so that one 
typically splits the whole time interval into subintervals of length r = At and 
then approximate 0{tH) acting on the initial condition at each step. Since 

^ra.A _( I akTH \ i,^B _( I 

V I )' V -bkTH I 

it makes sense to apply splitting methods of the form 

un+i = e^''"-^ e™"^ • • • e^^i^ e™^^ u^. (16) 

Observe that the evaluation of the exponentials of A and B requires only com- 
puting the products Hp and Hq^ and this can be done very efficiently with the 
FFT algorithm. 

Several methods with different orders have been constructed along these 
lines [131 [2T1 127] . In particular, the schemes presented in [13] use only m = r 
exponentials e^°' and e"^ ' to achieve order r for r = 4, 6, 8, 10 and 12. Fur- 
thermore, when the idea of processing is taken into account, it is possible 
to design families of symplectic splitting methods with large stability inter- 
vals and a high degree of accuracy [21 [3]. They have the general structure 
P{tH)K{tH)P~^{tH), where K (the kernel) is built as a composition (J16p 
and P (the processor) is taken as a polynomial. 



Although these methods are neither unitary nor unconditionally stable, they 
are symplectic and conjugate to unitary schemes. In consequence, neither the 
average error in energy nor the norm of the solution increase with time. In 
other words, the quantities ||m|P = u^u/N and u^ Hu/{2N) are both approx- 
imately preserved along the evolution, since the committed error is (as shown 
in Subsection 13.11 below) only local and does not propagate with time. The 
mechanism that takes place here is analogous to the propagation of the error in 
energy for symplectic integrators in classical mechanics [15]. In addition, the 
families of splitting methods considered here are designed to have large stability 
intervals and can be applied when no particular structure is required for the 
Hamiltonian matrix H. Furthermore, they can also be used in more general 
problems of the form q = Mip, p = —M2q, resulting, in particular, from the 
space discretization of Maxwell equations [3] . 

3 Analysis of symplectic splitting methods for time 
discretization 

In this section we proceed to characterize the family of splitting symplectic 
methods (116p . paying special attention to their stability properties. By inter- 
preting the numerical solution as the exact solution corresponding to a modified 
differential equation, it is possible to prove that the norm and energy of the 
original system are approximately preserved along evolution. We also provide 
rigorous estimates of the time discretization error that are uniformly valid as 
both the space and time discretizations get finer and finer. The analysis allows 
us to construct new methods with large stability domains such that the error 
introduced is comparable to the error coming from the space discretization. 

3.1 Theoretical analysis 

It is clear that the problem of finding appropriate compositions of the form (J16p 
for equation ()14p is equivalent to getting coefficients aj, bi in the matrix 

''^^''^=[-b^rH l)[o I )---[-b,rH l)[o I 

(17) 
such that K{tH) approximates the solution 0{tH), where 0{y) denotes the ro- 
tation matrix ()15|) . The matrix K{tH) that propagates the numerical solution 
of the splitting method (J17p can be written as 

^(^^) - V K,{tH) K,{tH) ) ' ^^^^ 

where the entries Ki{y) and Kiiy) (respectively, K2{y) and K^{y)) are even 
(repect., odd) polynomials in y G M, and detK{y) = Ki{y)Ki{y)—K2{y)K-i{y) = 
1. It is worth stressing here that by diagonalizing the matrix H with an ap- 
propriate linear change of variables, one may transform the system into N 
uncoupled harmonic oscillators with frequencies wi, . . . ^ijJn- Although in prac- 
tice one wants to avoid diagonalizing H^ numerically solving system (llSp by a 



splitting method is mathematically equivalent to applying the splitting method 
to each of such one-dimensional harmonic oscillators (and then rewriting the 
result in the original variables). Clearly, the numerical solution of each individ- 
ual harmonic oscillator is propagated by the 2x2 matrix K{y) with polynomial 
entries K.j{y) (j = 1,2,3,4) for y = tloj. We will refer to K{y) in the sequel 
as the propagation matrix, although other denominations have also been used 
[3l[23]. 

Moreover, for a given K{y) with polynomial entries, an algorithm has been 
proposed to factorize K{y) as (fTTj) and determine uniquely the coefficients aj, 
bi of the splitting method [Sj Proposition 2.3]. Thus, any splitting method is 
uniquely determined by its propagation matrix K{y). For this reason, in the 
analysis that follows we will be only concerned with such matrices K[y). 

When applying splitting methods to the system ()13p with time step size 
r, the numerical solution is propagated by {K{tH))'^ as an approximation to 
0{tH)^ = 0{nTH), which is bounded (with L^ norm equal to 1) indepen- 
dently of n. It then makes sense requiring that {K{tH))^ be also bounded 
independently of n > 1. This clearly holds if for each eigenvalue ujj of H, the 
corresponding 2x2 matrix K{y) with y = tujj is stable, i.e., if ||(i^(ra;j))"|| < C 
for some constant C > 0. 

In our analysis, use will be made of the stability polynomial, defined for a 
given K{y) by 

p{y) = \iiK{y) = \{Ki{y) + K^{y)). (19) 

The following proposition, whose proof can be found in [3], provides a charac- 
terization of the stability of K{y). 

Proposition 2 Let K{y) be a 2 x 2 matrix with det K{y) = 1, and p{y) = 
^tr K{y), with y € M. Then, the following statements are equivalent: 

(a) The matrix K{y) is stable. 

(b) The matrix K{y) is diagonalizable with eigenvalues of modulus one. 

(c) \p{y)\ < 1, and there exists a real matrix Q{y) such that 

Q{y)-^K{y)Q{y) = 0(<A(y)), (20) 

where 0{y) is the rotation matrix i f J5)) and (piy) = arccosp(y) € M. 

We define the stability threshold y* as the largest non negative real number 
such that K{y) is stable for all y € (— y*,y*). Thus, [K{tH))'^ will be bounded 
independently of n > 1 if all the eigenvalues of tH lie on the stability interval 
{—y*,y*), that is, if t p{H) < y^,, where p{H) is the spectral radius of the 
matrix H. For instance, if a Fourier-collocation approach based on N nodes is 
applied to discretize ([5]) in space, the spectral radius is of size p{H) = 0{N'^) 
(cf. eqs. ©-([H])), which shows that r must decrease proportionally to A^"^ as 
the number of nodes N increases. 



The stability threshold y* depends on the coefficients {oj, bi} of the method 
(J16p and verifies y* < 2m, since 2r?i is the optimal value for the stability thresh- 
old achieved by the concatenation of m steps of length r/m of the leapfrog 
scheme [7j. 

The stability of the matrix K{y) of a splitting method for a given y € M 
can alternatively be characterized as follows. 

Proposition 3 The matrix K{y) is stable for a given y gM. if and only if there 
exist real quantities (j){y),e{y),^{y), with 7(y) 7^ 0, such that 

cos(0(y)) + e{y) sin((/>(y)) -f{y) sin((/)(y)) 



7(y) 



■ sin(0(y)) cos((?;>(y)) - e(y) sin(</>(y)) 



Proof. If K{y) is of the form ([2T]) then, obviously, trii'(y) = 2cos((/)(y)) and 
thus cos((/>(y)) = p{y). Moreover, it is straightforward to check that (j20p holds 
with 



7(y)i/2 

-e(y)7(y)~^/^ 7(y)" 



Qiy) = ( ..::!i.-i/2 . .-1/2 ) , (22) 



so that K(y) is stable in that case. 

Let us assume now that K{y) is stable, so that from the third characteriza- 
tion given in Proposition [2l 0(y) = arccos(p(y)) G M, where p(y) is the stability 
polynomial. We now consider two cases: 

• p{y) = 1 (resp. p{y) = —1), so that K{y) (resp. —K{y)) is similar to 
the identity matrix, which implies that K{y) (resp. —K{y)) is also the 
identity matrix. In that case, (j2ip holds with e{y) = and 7(2/) = 1. 

• If p(y)^ 7^ 1, then sm{(j){y)) ^ 0, and we set 



2sin((/>(y)) sin((/>(y)) 

Since det(A'(y)) = 1, one has 

-K2{y)K3iy) = 1 - Ki(y)K4(y) = (1 + e{yf) sin(<A(y))^ 

which implies that 7(2/) 7^ and 

l + e(y)2 
K3{y) = - ,7 sin(0(y)). 

7(y) 



Notice that, for a given splitting method with a non-empty stability interval 
{—y*,y*), Proposition [3] determines two odd functions (j){y) and e(y) and an even 
function 7(2/) defined for y G (— y*,y*) which characterize the accuracy of the 
method when applied with step size r to a harmonic oscillator of frequency uj, 
with y = TLo. An accurate approximation will be obtained if |</'(y)— y|, |7(y) — 1|, 



and \e{y)\ are all small quantities. In particular, if the splitting method is of 
order r, then 

cPiy) = y + 0{f^^), e{y) = 0(/), ^{y) = 1 + 0{f) 

as 2/ — > 0. For instance, for the simple first order splitting e"^ e"^^ one has 

I y\ f I 0\ /l-y2y 



^(^) > 1 y V -y 1 y V -y i 

and one can easily check that 

(/)(y) = arccos(l - y^/2) = 2 arcsin(y/2) =y + 0{y^), 

e{y) = —L= = 0{y\ 
V4-y2 

7(y) = ^=L= = l + 0(y2). 
V4-y2 

It is worth stressing that (j2ip implies that (j20p holds with (5(y) given by 
(j22p . This feature, in particular, allows us to interpret the numerical result 
obtained by a splitting method of the form (|16p applied to (J14p in terms of the 
exact solution corresponding to a modified differential equation. Specifically, 
assume that qn + ipn = Un ^ u{tn) = exp{—itnH)uo is obtained (for i„ = rer, 
n > 1) as 



= KItHY 

Pn J \ PO 

If T p[H) < y*, then it holds that 

trivially verifies Un = exp(—in(p{TH))uo. In other words, Un is the exact solu- 
tion at tn = nr of the initial value problem 

i—u = Hu, u{0) = uq, (23) 

where H = -(J){tH) ~ H. With this backward error analysis interpretation at 
hand, it readily follows the preservation of both the discrete L^ norm of 

u = {-^{THy^l^ + i e{TH)-i{TH)~^l^) q + i -^{tH^/'^ p 

and the energy corresponding to (j23p . This implies that the discrete L^ norm 
of u = q + ip and the energy of the original system will be approximately 
preserved (that is, their variation will be uniformly bounded for all times tn)- 

3.2 Error estimates 

Our goal now is to obtain meaningful estimates of the time discretization error 
that are uniformly valid as N —^ oo and r — )• 0, that is, as both the space 
discretization and the time discretization get finer and finer. We know that, by 

10 



stability requirements, the time step used in the time integration by a sphtting 
method of system ([13]) must be chosen as r < y^/p{H), where the stabihty 
threshold must verify y^ < 2m for an m-stage splitting method. Since the 
Hermitian matrix H comes from the space discretization of an unbounded self- 
adjoint operator, the spectral radius p{H) will tend to infinity as A^ — ?• oo, and 
thus inevitably r — )• 0. It seems then reasonable to introduce the parameter 



e = Tp{H) 



(24) 



and analyze the time- integration error corresponding to a fixed value of ^ < y* . 
The fact that, for each y G {—y*,y*), (pOj) holds with Q{y) given by ([22]) 
implies that, for each n > 1, 



KiyY 



cos{n(j){y)) + e{y) sm{n(f){y)) 
l + e(y)2 



7(y) 



■ sm{n(j){y)) 



j{y)sm{n(f){y)) 
cos{n(p{y)) - e{y) s'm{ncl){y)) 



This will allow us to obtain rigorous estimates for the error of approximating 
uq by applying n steps of a splitting method with time-step r = t/n. 
Specifically, we have 



-itH 



K{tHY 



Po 



0{n4){TH)) 



m 

Po 



+ E{tH) 



sm{n(j){T H))qQ 
sm{n(j){T H))pq 



(25) 



where 0{y) denotes the rotation matrix ([15p and the 2x2 matrix E{y) is given 
by 

e(y) 7(y) - 1 



E{y) 



i{y) 



+ 1 



<y) 



(26) 



so that the following theorem can be stated. 



Theorem 4 Given uq = q^ + ipQ, let Un = q-n -\- iPn be the approximation to 



u{nT) 



-i nr H 



Uq obtained by applying n steps of length 



of a splitting 



method with stability threshold y^. Then one has 

\\un - u{nT)\\ < {np{9) + 7/(6*)) ||tio| 
(in the Euclidean norm), where 



sup |0(y) -y\, 

0<y<9 



Proof. From ( 1251 ). we can write 

\\Un — u{nT 



HO) 



sup ||£' 
0<y<8 



Qn - q{tn) 
Pn - P(tn) 



< 



{0{n(P) - 0{nTH)) 



qo 

Po 



+ \\E(tH) 



sm{n<j))qo \ 
sm{n(l))po J 
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where, for clarity, (j) = 4'i'^H). For the first contribution we have 



{0{n(t)) - 0{nTH)) 






(e-^"'^-e-^"^^)uo| 



|g-mT_ff^ _ ^-in(4)-TH)\ 



< iia-e- 



-in{<l)—TH)\ 



ml 



since H is Hermitian. Now 



-in{(l)—TH)s 



in I e 



< n I \\e-M<P-rH)s 
Jo 

= n max \(h(TUJi) — TUJi\ < n n(9) 
l<j<N '^^ J' ■" 



TH)ds\\ 

TH\\ds = nU-TH\\ 



As for the second contribution, one has 



sm{n(f>)qo 
sin(n0)po 



{n(j))uQ\\ < ||sin(n<^)|| ||uo|| < ||tio| 



whereas 



\E(tH) 



sm ni 



max WEiTiOn] 



< i^ie). 



and thus the proof is complete. ■ 

Notice that the error estimate in the previous theorem does not guarantee 
that, for a given t, the error in approximating e~* uq by applying n steps of 
the method is bounded as p{H) — )■ oo. As a matter of fact, since r = t/n must 
satisfy the stability restriction 9 = Tp{H) < y^,, so that n > t p{H)/y^, one has 
that n (and hence the error bound above) goes to infinity as p{H) -^ oo. This 
can be avoided by estimating the error in terms of ||i7Mo|| in addition to |Iuo||. 
The assumption that ||ii^tio|| can be bounded uniformly as the space discretiza- 
tion parameter A'^ — >■ oo, implies that the initial state ■ip{x,0) of the continuous 
time dependent Schrodinger equation is such that d'^'il){x, 0) is square-integrable. 
The converse will also be true for reasonable space semi-discretizations and a 
sufficiently smooth potential V{x). 

More generally, the assumption that '(/'(x, 0) has sufficiently high spatial 
regularity (together with suitable conditions on the potential V{x)) is related 
to the existence of bounds of the form ||i7'^no|| < C^ that hold uniformly as 
p{H) -^ oo. In this sense, it is useful to introduce the following notation: 

• Given /c > 0, we denote for each u € C^ 



\u\\k ■■-- 



\H^u\ 



a 771-stage splitting method with stability threshold y*, given /c > 
and € [0,2/,,) we denote 



For 



Pk{0) = sup 
o<y<e 



ct>{y) 



{0/y)\ 



MO) = sup ||i^(y)|l(^/y)'=. 
o<y<e 



(27) 
(28) 



uo\ 
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Clearly, /ifc(0) and z^fc(0) are bounded if and only if the method is of order 
r > k. 

We are now ready to state the main result of this section. 

Theorem 5 Given uq = qq + ipo and t £ M'^", let n he such that r = t/n = 
9/p{H) (with 6 < y^), and let Un = Qn + ipn be the approximation to u{t) = 
g-jtH^^ Qlldined by applying n steps of length t of a r-th order splitting method 
with stability threshold y*. Then, for each k € [0, r], 



tfikiO) ll^to llfc+i + i^k{0)\\uo\\k 
p{HY 



\un-u{t)\\< ^j—j: . (29) 



Proof. We proceed as in the proof of Theorem HI First we bound 

||e-^""^wo - e-^"'^("^)no|| < T^+^ \\{THy^^\l - e-*"W-^)-^))|| ||uo||fc+i 

te^ 



- -^(J^\\i^H)-^^\<t^{TH)-TH))\\\\u,U+i. 



with 



\\{TH)-'^-'{<t,{TH)-TH)\\ = ^max^|(ra;,)-'''-'(0(r^,) -to;,)! < ^^k[e)/e\ 
Then the second term in (1251) verifies 



r^ \\{THy^E{rH)s\n{n(t>{rH))H^u4 < -^\\{THr''E{T H)\\ ||«olU, 

p{Hf 

and 

\\{tH)-'E{tH)\\ = max ||(™,)-'^(™.)ll < Me)/e\ 

1<J<N 

from which ()29p is readily obtained. ■ 

Some remarks are in order at this point: 

1. Recall that the estimate in Theorem [1] shows the behavior of the space 
discretization error (of a spectral collocation method applied to the ID 
Schrodinger equation) as the number A^ of collocation points goes to infin- 
ity. Our estimate (j29p shows in turn the behavior of the time discretization 
error as A^ ^' co, provided that r = 9/p{H) with a fixed 9 < y* . In that 
case, it can be shown that p{H)~^ < L N~'^ uniformly for all N, and thus 
the error of the full discretization admits the estimate 

^ (C(l + t) max ||9f +2V;(-,0|| + i(t/^fcWll^olU+i + ^fcWIIuolU)) • 

Notice the similarity of both the space and time discretization errors 
(||no||fc+i = ll-ff^'^^uoll is a discrete version of a continuous norm ||V'(-,0)||fc+i 
which is equivalent to the Sobolev norm ||5^'^"''^V'(")0)II)- 
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2. Given a splitting method with stabihty threshold y* of order r for the 
harmonic oscillator, consider firiO) and Ur{9) in ([27j) - (p8|) for a fixed 
6 < y*. If instead of analyzing the behavior as A^ increases of the time dis- 
cretization error committed by the splitting method when applied with 
T = 6/p{H), one is interested in analyzing the error with fixed H and 
decreasing t < 6/p(H), one proceeds as follows. Since by definition 

Pr{rp{H)) < pM i^^] ' , Mrp{H)) < MO) ^ ^^^^^ 



6 
then reasoning as in the proof of Theorem [5l one gets the estimate 

II /^Mi / */^r(6') ||uo||r+l + t'r(6')||uo||r r 
\\Un - U[t) II < T . 



From a practical point of view, (|29p can be used to obtain a priori error 
estimates just by replacing the exact p{H) by an approximation (obtained for 
instance with some generalization of the power method), or by an estimation 
based on the knowledge of bounds of the potential and the eigenvalues of the 
discretized Laplacian. 

The error estimates in Theorem [5] provide us appropriate criteria to con- 
struct splitting methods to be applied for the time integration of systems of the 
form (I13p that result from the spatial semi-discretization of the time dependent 
Schrodinger equation. Such error estimates suggest in particular that different 
splitting methods should be used depending on the smoothness of initial state 
in the original equation. Also, Theorem [5] indicates that for sufficiently long 
time integrations, the actual error will be dominated by the phase errors, that 
is, the errors corresponding to pkiO)- 

4 On the construction of new symplectic splitting 
methods 

Observe that when comparing the error estimates in Theorem [5] for a given 
A; > corresponding to two methods with different number of stages m and m! 
respectively, one should consider time steps r and r' that are proportional to 
m and m! respectively. In this way, the same computational effort is needed for 
both methods to obtain a numerical approximation of u{t) for a given i > 0. It 
makes sense, then, to consider a scaled time step of the application with time 
step r of a m-stage splitting method to the system ()13p . This can be defined 
as 

»' . 1 = I^, (30) 

m m 

so that the relevant error coefficients associated to the error estimates in The- 
orem [5] are pk{S'm) and Vk{9'm). 

The task of constructing a splitting method in this family can be thus pre- 
cisely formulated as follows. 
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Problem. Given a fixed number m of stages in ( (J7| j, and for prescribed values 
of k > and scaled time step 9' £ (0,2), design some splitting method having 
order r > k and stability threshold y^ > 6'm, which tries to optimize the main 
error coefficient fikid' in) while keeping i'k{9'm) reasonably small. 

We have observed, however, that trying to construct such optimized meth- 
ods in terms of the coefficients of the polynomial entries Kj{y) [j = 1, 2, 3, 4) of 
the propagation matrix K[y) leads us to very ill-conditioned systems of alge- 
braic equations. That difficulty can be partly overcome by taking into account 
the following observations: 

• The functions \4>{y)/y — 1| and ||i?(y)|| (y G {—y*-,V*)) determining the 
error estimates in Theorem [5] uniquely depend on two polynomials: the 
stability polynomial p{y) given in (I19p and 

Indeed, from one hand, 4>{y) = arccos(p(y)), so that \4>{y)/y — 1| uniquely 
depends on p{y). On the other hand, according to Proposition [3l 

<l{y) = (1 + T,Ky)) sin((A(y)), where 5{y) = ( 7(2/) + 



and one can get by straigthforward algebra that ||S(y)|| is (in Euclidean 
norm) 



\E 



\ 



,„) ,,M,,L,),M! 



(and thus ||i^(y)|| = ,/5{y) + 0{5{y)) as 5{y) ^ 0). 

Given an even polynomial p{y) and an odd polynomial g(y), there exist a 
finite number of propagation matrices K{y) such that (fT9]) and (f3T]) hold. 
Indeed, the entries of such stability matrices are of the form 

Ki{y) = p{y) + d{y), K2{y) = q{y) + e(y), 
K^{y) = -q{y) + e(y), K^iy) = p{y) - d{y), 

where d{y) and e(y) are respectively even and odd polynomials satisfying 

p{y?+qiyf-l = d{y)^ + eiyf. (32) 

It is not difficult to see that there is a finite number of choices for such 
polynomials d{y) and e(y). The ill-conditioning mentioned before seems 
to come mainly from the ill-conditioning of the problem of determining 
d{y) and e(y) from prescribed polynomials p{y) and q{y). Obviously, a 
necessary condition for the existence of such polynomials d{y) and e(y) 
with real coefficients is that p(y)^ -|- q{y)'^ — 1 > for all y. It is also 
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straightforward to see that, for a rth order method, d{y) = 0{y^~^^) and 
e{y) = C'(y''+^) (as y -^ 0), and thus 

P(y)' + 9(y)'-i = 0(y'^+'). (33) 

In addition, if the method has stabihty threshold y^, > 0, then there exists 
< yi < ■ ■ ■ < yi < y^ such that (p{yj) = jtt, and thus 

p{yj) = {-iy, p'(y,) = 0, q{yj) = 0, for j = l,...,L (34) 

For simphcity, we restrict ourselves to the construction of m-stage methods 
of even order r that are intended to have small values of /ir(0'r7i) and Vr{6'nn) 
for a prescribed scaled time step 9' E (0,2). When designing such a method, 
we follow several steps: 

1. First find two polynomials p{y) and q{y) with small value of iir{0'm) + 
\vr{0'm) (for some A < 1) among those satisfying the following three 
conditions: 

(a) There exist yj ~ jvr [j = 1, . . . ,/) with In < O'm < (/ + l)7r such 
that dal) holds; 



(b) p{y) = cos(y) + 0(y^+i), q{y) = sm{y) + 0{y'+^), and ([MD as y ^ 0; 

(c) p{yf + q{yf - 1 > for ah y G E. 

2. Find all possible pairs {d{y),e{y)) of real (even and odd respectively) 
polynomials satisfying ([32]) . and for each pair {d{y),e{y)), construct the 
corresponding 2x2 matrix K{y). 

3. Apply the algorithm given in [3] to each of the matrices K{y) obtained in 
the previous step. In this way we will get the vector of coefficients {ai,bi) 
of all the splitting methods having a progagation matrix K{y) satisfying 
(fT9]) and ([3T]) for the pair of polynomials {p{y),q{y)) determined in the 
first step. Since Theorem [5] gives exactly the same error estimate ()29p for 
all the splitting schemes obtained in that way, we choose (with the aim 
of reducing the effect of round-off errors) one that minimizes 

m 



This procedure has been applied to construct several splitting methods of 
different orders r, number of stages m and scaled time steps 9' . We collect in 
Table [T]the relevant parameters of some of them, whereas the actual coefficients 
aj,hj can be found at www. gicas .uji .es/Research/splittingl .html. As a 
matter of fact, all the methods have m -\- 1 pairs of coefficients a^, 6j, but 
hm+i = 0. In consequence, the last stage at a given step can be concatenated 
with the first one at the next step, so that the overall number of stages is 
m. This property is called FSAL (first-same-as-last) in the numerical analysis 
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m 


r 


e' 


y*/m 


E,(l«.l + I^.l) 


flri^'m) 


Ur{e'm) 


10 


6 


1 


1.1617 


4.022 


0.0009341 


0.0372 


20 


16 


1 


1.0456 


3.0553 


0.000611028 


0.0258433 


30 


24 


1 


1.0246 


3.19658 


0.0000841871 


0.0373544 


30 


6 


1.4 


1.41876 


3.0921 


0.0000518519 


0.0131295 


30 





1 


1.1411 


3.04948 


2.91902 • 10-13 


2.28673 • 10-9 


30 





0.75 


1.027 


3.44381 


1.2545 • 10^1^ 


5.96706 • 10-1^ 


30 





0.5 


0.937874 


3.84442 


7.96031 • 10^24 


6.66693 • 10-is 


40 





1 


1.15953 


3.21986 


1.06301 • 10-15 


1.07587 • 10-12 



Table 1: Relevant parameters of several new splitting methods of order r es- 
pecially designed to integrate with scaled time step 9' = Tp{H)/m the semi- 
discretized Schrodinger equation. Here p{H) is the spectral radius of the matrix 
H, m is the number of stages, y* stands for the stability threshold and Pri^), 
Vr{9) are the coefficients appearing in the error estimate (I29p . 



literature. According to the previous comments, the new schemes are aimed at 
integrating equation (J13p under very different conditions of regularity. 

The first three methods in Table [U are designed to be applied with the same 
scaled time-step 9' = Tp{H)/m = 1, and thus the three of them have the same 
computational cost. The order r of the methods is increased by adding more 
stages, while keeping reasonably small error coefficients iJir{9) and Vr{9)- This 
will be advantageous, according to Theorem [5l for sufficiently regular initial 
states. Alternatively, one may want to use the additional number of stages to 
reduce the computational cost while keeping the same order r = 6. This can 
be illustrated with the fourth method in Table [H which has been optimized 
for scaled time-step 9' = 1.4, and thus is substantially cheaper than the first 
method (optimized for 9' = 1), while having smaller error coefficients Pq{9) and 

We now turn our attention to the methods of order zero in Table [H which 
according to Theorem [5l are the methods of choice for very low regularity 
conditions. Although they have comparatively smaller error coefficients than 
the methods of order r > 1, one should bear in mind that the error estimates 
(I29p for r >k > 1 decrease with p{H)~ as the spectral radius p{H) increases, 
while for methods of order r = the same estimate holds independently of the 
size of p{H). Comparing the first three methods of order r = and in = 30, 
we see that, not surprisingly, the accuracy of the methods can be improved by 
increasing the computational cost (by considering methods optimized for lower 
values of 9') . This is analogous to increasing the accuracy of the application of a 
Chebyshev polynomial of degree m = 30 by decreasing the time-step size. Now, 
by comparing the first 30-stage method of order with the method with m = 40 
and order 0, we see that the accuracy can be increased also by increasing the 
number of stages from tti = 30 to ?ti = 40 while keeping the same computational 
cost (with 9' = 1). This is similar to increasing the accuracy of Chebyshev 
approximations, while keeping the same computational cost, by increasing the 



17 





t \ojj\nj + Uj, 








(myj) 


, Uj = \\E{myj)\\ with 


\y^\ 


TUJj 

m 


P{H) 



degree of the polynomial from m = 30 to jtt, = 40. 

Recall that, if wi, . . . ,ujn are the eigenvalues of H , numerically integrating 
p3|) by a splitting method is mathematically equivalent to applying the splitting 
method to N uncoupled harmonic oscillators with frequencies ujj. Particular- 
izing the proof of Theorem O to this case, it is quite straightforward to con- 
clude that, when integrating the system with scaled time step 6' (that is, with 
r = mO'/p{H), where m is the number of stages of the scheme,) the relative 
error made in each oscillator can be bounded by 



where 

and thus \yj\ < 9'. When such a system of harmonic oscillators originates from 
a continuous problem possessing a high degree of regularity, the highest fre- 
quency oscillators have much smaller amplitude and thus can be approximated 
less accurately than the lower frequency oscillators without compromising the 
overall precision. For lower regularity conditions, the overall precision will be 
more affected by the accuracy of the approximations corresponding to higher 
frequency oscillators. 

With the aim of illustrating the relative error made in each harmonic oscilla- 
tor, in Figure[T]we represent (in double logarithmic scale) \(l){my) / {my) — 1\ and 
||-E'(r7T,y)|| (which are even functions of y) for y £ [0, 9'] for some of the methods 
collected in Table [U identified by appropriate labels indicating their respective 
number of stages, order and scaled time step {m, r, 9'). Observe that both func- 
tions of y exhibit a similar behavior for each splitting method, although the 
values taken by the second one are several orders of magnitude smaller, since 
the methods are designed to minimize mainly the phase error coefficient ^k{9)- 
The order r of each of the methods is reflected in the slope of the curves as y 
approaches 0. 

We also include in Figure[T]the 12th order 12-stage scheme presented in jlSj . 
which is perhaps the most efficient when applied to harmonic oscillators among 
those (non-processed) splitting methods currently found in the literature. We 
denote it by GM12. It has a relative stability threshold y^/12 = 0.2618, so that, 
strictly speaking, it should be used with 9' = Tp{H)/12 < 0.2618 to guarantee 
stability. However, it seems in practice that the method can be safely used with 
9' = 0.932 (for a larger value of the scaled time step 9', the method becomes 
very unstable), because \\p{y)\ - 1| < lO"*' provided that \y\/12 < 0.932183. 
The theoretical instability for 9' G (0.2618, 0.932183) is only relevant after a 
very large number of steps, and reveals itself as resonance peaks (which are 
clearly visible in the graph of ||ii^(my)|[ in Figure [Ufor A; = 2, 3) near the values 
9' = k7r/12, k = 1,2,3. 

We can see in Figure [T] that GM12 is less accurate than the 30-stage 24th 
order method for the whole frequency range, and thus the former will show a 
poorer performance than the later for any regularity conditions. If the ampli- 
tudes at higher frequencies decrease fast enough, the 24th order method will 
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0.500 1.000 



Figure 1: Graphs of \4>(my) / (my) — 1| (top) and ||i?(?TT,y)|| (bottom) for some of 
the TTi-stage sphtting methods collected in Table [1] and the 12th-order scheme 
GM12. Each new splitting method is identified by the triad {m, r, 6'), indicating 
its number of stages m, order r and scaled time step 6', as in Table [TJ The 
resonances associated with the instability of GM12 at y ~ fcvr (A; = 2, 3) are 
visible, especially in the second graph. 



give very accurate approximations of u{t) = e~**^ uq at a relatively low cost 
(since 6' = 1). The 6th order method with m, = 30 stages gives correct approx- 
imations for all harmonic oscillators within the range y € [—1.4, 1.4], and hence 
it is expected to give excellent approximations under mild regularity with a 
comparatively lower cost than the 24th order method {9' = 1.4 compared with 
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9' = 1). Clearly, the methods of order will be the right choices for low reg- 
ularity conditions, since the corresponding phase errors \4>{my) / {my) — 1| are 
uniformly bounded for all y £ [—1, 1]. Among them, the method with m = 30 
and 6' = 1 can be accurate enough in many practical computations. If more 
precision is required, one can either consider the method with 6' = 0.75 (with 
result in a 25% increase of the computational cost), or use the method with 
771 = 40 and 9' = 1, without any increase of the computational cost (at the 
expense of having less frequent output). 

5 Numerical examples 

The purpose of this section is twofold. On the one hand, since the symplectic 
splitting methods we have presented here to approximate e~**^«o involve only 
products of the matrix H with vectors, it makes sense to compare them with 
other well established schemes of this kind, such as the Chebyshev and Lanczos 
methods. Although a thorough comparison with the family of splitting methods 
proposed in this work will be the subject of a forthcoming paper [4], we include 
here some results which show that the new schemes are indeed competitive for 
evaluating exp{—itH)uQ, at least in the example considered. 

On the other hand, since Theorem [5] provides a rigorous a priori estimate 
on the error committed when using a splitting method of the form (I17p in the 
time integration of equation ([9]), it is interesting to check how this theoretical 
error estimate behave in practice for some of the methods constructed here. 

5.1 A preliminary comparison with Chebyshev and Lanczos 
methods 

As is well known, Chebyshev and Lanczos methods provide high order poly- 
nomial approximations to e~ uq requiring only matrix- vector products. The 
former is neither unitary nor symplectic, whereas the later is unitary, but sym- 
plectic only in the Krylov subspace (which changes from one time step to the 
next). 

To carry out this comparison we choose the very simple example previ- 
ously considered in [22]. The problem consists in approximating y = e~ u, 
where u is a random vector of unit norm and A is the tridiagonal matrix 
A = ^ tridiag(— 1, 2, — 1) of dimension 10000. The eigenvalues of A are con- 
tained in the interval [0, 2a;]. 

We have implemented the Chebyshev and Lanczos algorithms in the usual 
way (see [22]) with the particularity that, since the range of values for the 
eigenvalues is known, both the Chebyshev and the new splitting methods are 
used with a shift to the midpoint of the spectrum. In other words, y = 
^-luii ^-i{A-uii)y^ with / the identity matrix. This shift allows us to take 
p{A — ujI) ~ oj. 

Figure [2] shows the error, ||y — yapll for different degrees m of the polynomials 
used and for uj = 15, 20, 30, 40. Here y is computed numerically to high accuracy 
and yap corresponds to the approximate solution obtained by each scheme. 
Each particular value of m in the Lanczos and Chebyshev methods corresponds 
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Figure 2: Error, \\y — yap\\, versus degree of the polynomial, m, for approxi- 
mations to y = e~ V where v is chosen as a random vector of unit norm and 
A is the tridiagonal matrix A = ^ tridiag(— 1,2, —1) of dimension 10000. Here 
y is computed numerically to high accuracy and yap correspond to the approx- 
imate solutions obtained by each method. Results corresponding to Lanczos 
(crosses), Chebyshev (small circles) and several new splitting {m,r,6') methods 
(big circles) are depicted. 

to a different mth-order polynomial approximation (denoted by small crosses 
and circles, respectively). We clearly observe that, for this irregular problem, 
the Lanczos method converges to the optimal Chebyshev method, the main 
difference between both schemes being the number of vectors to be kept in 
memory. 

To apply the new splitting methods, we notice that for this problem the 
time step r = 1, and the corresponding spectral radius p ^ oj. Therefore we 
shall consider splitting methods whose value of 6' given by (f30ll satisfies 

m m 
In other words, for each co the method {m,r,9') is such that mO' >uj. 
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From the graphs of Figure [21 it is clear that, for each value of w, we can al- 
ways select one particular splitting method (big dots) which outperforms both 
Chebyshev and Lanczos. High-order splitting methods show a worst perfor- 
mance than schemes of order zero for this problem, since they require typically 
a higher degree of regularity. 

5.2 The Poschl-Teller potential 

We next illustrate the error estimate provided by Theorem for the class of 
splitting methods proposed here. We also compare the error in the time inte- 
gration with the error coming from the space discretization for different values 
of the mesh size N. For that purpose we choose a well known anharmonic 
quantum potential leading to analytical solutions and consider, for clarity, only 
the 30-stage splitting methods of order six and order zero collected in Table [TJ 
Specifically, we consider the Poschl-Teller potential 

a^ A(A - 1) 



V{x) 



2/i cosh (ax) 



with A > 1. It has been frequently used in polyatomic molecular simulation and 
is also of interest in super symmetry, group symmetry, the study of solitons, etc. 
[H [ini EO]- The parameter A gives the depth of the well, whereas a is related 
to the range of the potential. The energies are 



Q,2 

En = - — (A - 1 - nf, with < n < A - 1. 



We take the following values for the parameters (in atomic units, a.u.): 
the reduced mass n = 1745 a.u., q = 2, A = 24.5 (leading to 24 bounded 
states), X G [—5, 5], and assume the system is periodic. The periodic potential is 
continuous and very close to differentiable. For A^ = 128 we have p{H) ~ 0.635, 
whereas for N = 256 one gets p{H) ~ 1.85. 

We take as initial condition the Gaussian function, Tp{x, 0) = a e~^ ^ , where 
(J is a normalizing constant. With 6 = 3 the function and all its derivatives of 
practical interest vanish up to round off accuracy at the boundaries. The initial 
conditions contain part of the continuous spectrum, but this fact does not cause 
any trouble due to the smoothness of the periodic potential and wave function. 
As an illustration, some of the corresponding values of ||?io||fc for A^ = 128 
are: ||no||i = 0.629909, ||no||6 = 0.0722513, ||no||7 = 0.0478388. They decrease 
only moderately for the first values of k (before they increase again due to 
the contributions coming from higher energies). The corresponding values for 
A^ = 256 are quite similar to the previous ones. For this problem, both large 
and very small spatial errors are expected from spectral methods, depending on 
the mesh employed. It is then useful to have different methods with large values 
of 9' when low accuracy is desired and smaller values of 9' for high accuracy. 

We integrate for t € [0, 128r] with T = 333 and measure the 2-norm er- 
ror in the discrete wave function, ||nex(2*r) — Uap(2*r)||, for i = 0, 1, . . . , 7. 
The values Uex ai'e computed using the same spatial discretization and an ac- 
curate time integration (using a very small time step), whereas Uap stand for 
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Figure 3: Error in time integration (solid lines), error bounds from (j29p (dotted 
lines) for methods (30, 6, 1.4) and (30, 0, 1), and the spatial error (dashed lines), 
along the interval t G [0, 128T] for the Poschl-Teller potential. We have also 
included in the right panel the results obtained with the 30-stage method of 
order zero and 9' = 0.75, (30,0,0.75). 



the numerical approximations obtained with splitting methods. For a given 
spatial discretization, we choose the time step for each of the new methods 
such that T < mO' /p{H). In consequence, a period T has to be divided into 
M steps such that M = T/t > Tp[H)/{m9'). In particular, for the 6th-order 
method (30,6,1.4) and N = 256, since p{H) ~ 1.85, we take M = 15 > 
(333 X 1.85)/(30 x 1.4) (each period T requires 15 steps, 450 stages or 1800 
FFTs calls). 

The results obtained are shown in Figure [3l Solid lines represent the error 
with respect to the exact solution for the same spatial discretization, whereas 
dashed lines correspond to the total error with respect to the exact solution 
obtained with a finer mesh (it is the sum of the spatial error and the error from 
the time integration). Dotted lines are obtained with the estimate (j29p . 

We observe that the spatial error decreases exponentially with N due to 
the smoothness and periodicity of the problem. To estimate this spatial error 
we take the results obtained with A^ = 512 and an accurate time integration 
as the exact solution, and compare with the solution computed up to a high 
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accuracy for A^ = 128 and N = 256. For N = 128 the spatial error dominates 
the total error, so that the most convenient time integration scheme is one 
able to provide such accuracy with a large time step. These requirements are 
fulfilled by the (30, 6, 1.4) method, especially designed to be used with 9' = 1.4, 
whereas scheme (30, 0, 1) gives us higher accuracy than necessary and with more 
computational cost. Method (30,6, 1.4) can be used with a time step r about 
a 40% larger than method (30, 0, 1), and thus its computational cost is reduced 
approximately by this factor . 

However, for N = 256 the spatial error reaches nearly round off accuracy, 
and it could be convenient to employ methods able to provide this accuracy with 
the minimal computational cost. Notice that in this case the error committed by 
the 30-stage method with 0' = 1, (30, 0, 1), is still larger than the spatial error. 
In consequence, it makes sense integrating in time with a method specially 
designed to be used with a smaller time step. Thus, in particular, we reach 
round off accuracy with the 30-stage method (30, 0, 0.75) {6' = 0.75) which is 
nearly twice more expensive than the method with 6' = 1.4. 

We have also performed here the time integration with the 12th-order scheme 
GM12, which in the case of A^ = 128 requires a scaled time step of 6' = 0.49 
to give a precision similar to that obtained by (30,6, 1.4) with 9' = 1.4. With 
A^ = 256, GM12 must be applied with 6' = 0.19 to achieve the precision ob- 
tained by (30,0,0.75) with 6' = 0.75, thus requiring approximately four times 
more FFT calls. 

6 Concluding remarks 

The time integration of the Schrodinger equation previously discretized in space 
has been extensively studied in the literature. This is essentially equivalent to 
approximate u{t) = e~**^u(0), where H is a real symmetric matrix and u(0) 
represents the discrete wave function. In this work we propose using sym- 
plectic splitting integration methods to get this approximation. The main 
difference with standard polynomial approximations is that in the products 
Hv = HIie{v) + iHlin{v), the real and imaginary parts are computed sequen- 
tially instead of simultaneously (i.e., the computation of the real part is used 
in the computation of the imaginary part and vice versa in consecutive stages). 
These schemes are conjugate to unitary methods, so that the errors in norm 
and energy do not grow secularly [3]. 

To carry out the integration, one divides the whole time interval into n steps 
of length T = t/n and applies an m-stage method at each time step. The total 
computational cost of the method is measured by the product n m instead of 
m. The analysis carried out in this paper allows us, in particular, to construct 
a particular symplectic splitting scheme of the form (J16p which minimizes the 
total cost nm, given a a prescribed tolerance, the spectral radius p{H) of the 
corresponding Hamiltonian matrix H and the norm of its action on the initial 
condition, ||ii^'^u(0)||. We have observed that the optimal methods in this sense 
have relatively large values of m. We can choose the most appropriate method 
for each problem, i.e. the method, {m,r,9'), with the largest value of 6' which 
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provides the desired accuracy for a given problem. 

The error analysis of splitting methods provided here allows one to get a 
priori bounds on the propagating error when numerically integrating with a 
given time step which are comparable to similar estimates for the space dis- 
cretization error. Moreover, it permits to construct new classes of schemes with 
a large stability interval specifically designed to be used with a certain (large) 
time step in such a way that the accuracy is similar to the spatial discretization 
error for a given space regularity. The main ingredients in the process are again 
the values of p{H), \\H''u()\\, and the estimate provided by Theorem [H The 
numerical examples considered illustrate the validity of our approach. In par- 
ticular, they show that there are methods in this family which are competitive 
with other standard procedures, such as Chebyshev and Lanczos methods. 

By following this procedure it is indeed possible to generate a list of inte- 
gration schemes specifically designed to be used under different regularity con- 
ditions on the initial state and the Hamiltonian matrix which involve in each 
case an error comparable to that coming from the spatial discretization. It is 
our purpose in a forthcoming paper [5] to elaborate an algorithm in such a way 
that, given a prescribed tolerance, an initial state uq and a Hamiltonian ma- 
trix H, automatically selects the most efficient time integration method in this 
family fulfilling the requirements supplied by the user. Moreover, we will also 
carry out a detailed numerical study of this family of splitting methods and the 
proposed automatic algorithm in comparison with the Chebyshev polynomial 
expansion scheme and the Lanczos iteration method. 
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